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We present measurements of the bispectrum of dark matter halos in numerical simulations with 
non-Gaussian initial conditions of the local type. We show, in the first place, that the overall effect 
of primordial non-Gaussianity on the halo bispectrum is larger than on the halo power spectrum 
when all measurable configurations are taken into account. We then compare our measurements 
with a tree-level perturbative prediction finding good agreement at large scale when the constant 
Gaussian bias parameter, both linear and quadratic, and their constant non-Gaussian corrections are 
fitted for. The best-fit values of the Gaussian bias factors and their non-Gaussian, scale-independent 
corrections are in qualitative agreement with the peak-background split expectations. In particular, 
wo show that the effect of non-Gaussian initial conditions on squeezed configurations is fairly large 
(up to 30% for Jnl = 100 at redshift z = 0.5) and results from contributions of similar amplitude 
induced by the initial matter bispectrum, scale-dependent bias corrections as well as from nonlinear 
matter bispectrum corrections. We show, in addition, that effects at second order in /nl are 
irrelevant for the range of values allowed by CMB and galaxy power spectrum measurements, at 
least on the scales probed by our simulations. Finally, we present a Fisher matrix analysis to assess 
the possibility of constraining primordial non-Gaussianity with future measurements of the galaxy 
bispectrum. We find that a survey with a volume of about 10 h~^ Gpc^ at mean redshift z ~ 1 
could provide an error on /atl of the order of a few. This shows the relevance of a joint analysis of 
galaxy power spectrum and bispectrum in future redshift surveys. 
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I. INTRODUCTION 

Recent years have witnessed an intense activity related to the possibility, offered by future observations, 
to detect a non-Gaussian component in the primordial density perturbations. In fact, the detection 
of a non-vanishing primordial bispectrum would have profound implications for our understanding of 
the inflationary mechanisms, possibly nding out the simplest model of canonical, single-field, slow roll 
inflation [1]. 

On the theoretical side, primordial non-Gaussianity (PNG) has been recognized as a relevant "predic- 
tion" of early Universe models which has the potential to shed light on inflaton interactions or multiple 
fields scenarios (see [2-4] for recent reviews and [5, 6] for an effective field theory approach). On the 
observational side, the Planck satellite [7] will soon improve significantly the constraints on non-Gaussian 
parameters from the Cosmic Microwave Background (CMB) bispectrum measurements currently provided 
by the WMAP mission [8]. 

While the CMB clearly represents the most direct window on the primeval perturbations, the effect of 
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non-Gaussian initial conditions on the large scale bias of halos and galaxies discovered by [9] has turned 
the galaxy power spectrum into a powerful probe of PNG of the local kind [10], or, more gcnerically 
but to a lesser degree, of an initial bispcctrum assuming large values in the limit of squeezed triangular 
configurations such as in the folded model [11-13] or quasi-single field models [14, 15]. This discovery has 
been followed by a number of studies investigating this effect in details [16-34] and has motivated several 
different groups to runs large suites of N-body simulations seeded with non-Gaussian initial conditions 
[21, 24, 29, 30, 32, 35-37]. Remarkably, observations of the large scale structure can now place limits 
on the local non-Gaussian parameter /jvl similar to those from the CMB [17, 38, 39], while even better 
constraints are expected from upcoming large- volume galaxy redshift surveys [40-45], especially if the 
stochasticity induced by sampling variance and shot noise can be suppressed [46-49]. 

The interest in large-scale structure observations as a probe of the initial conditions has, however, a 
much longer history. In fact, it has long been expected that the main effect of PNG on the large scale 
mass and galaxy distribution consists in an additional contribution to the matter bispectrum, linearly 
evolved to present time by gravitational instability (see [50] and references therein). Focusing on this 
effect, [51, 52] have shown that future galaxy surveys will be able to provide constraints on the non- 
Gaussian parameters which surpass the best CMB limits. These results follow from the fact that future 
three-dimensional redshift surveys will provide a larger number of observable modes than two-dimensional 
CMB observations. Moreover, the primordial contribution to the matter bispectrum is to the large scale 
structure what the CMB bispectrum is to the temperature anisotropics. Hence, the impact of non- 
Gaussian initial conditions of the local type on the galaxy power spectrum should not overshadow the 
sensitivity of the galaxy bispectrum to PNG. In fact, the latter is much more sensitive to the amplitude 
and configuration shape of the three-point function of primordial curvature perturbations. Therefore, 
while the galaxy power spectrum might be competitive with the CMB bispectrum as far as local non- 
Gaussianity is concerned, this does not happen for other models of non-Gaussianity. In more general 
terms, a complete study of large-scale structure data will naturally involve the combined analysis of 
power spectrum and bispectrum measurements (see, for example, [53, 54] in the context of cosmological 
parameter determination). This work represents a first step in the direction of such an analysis for 
non-Gaussian initial conditions. 

In this work, we focus on the local model of primordial non-Gaussianity [55-57] . In this specific case we 
expect, for the halo bispectrum the combination of effects due to scale-dependent corrections to the both 
linear and quadratic halo bias and the effect of the primordial bispectrum component. First attempts at 
incorporating large-scale bias corrections in a description of the halo or galaxy three-point function can 
be found in [58, 59]. These works relied however on a local model for galaxy bias or high-peak statistics 
leading to equivalent results. It is now clear that such a prescription does not correctly describe two-point 
statistics (see, for instance, [32, 33]). A simple description that takes into account nonlocal corrections 
in terms of a multivariate bias expansion is proposed by [25], developing the earlier results of [17-21] 
to include nonlinear bias contributions. They also provide a partial expression for the halo bispectrum 
including a few relevant terms. A complete tree- level expression for the halo bispectrum based on the 
multivariate halo bias expansion is derived instead in [60] . This is essentially the model that we compare 
to numerical simulations in this work. 

Previous measurements of the halo bispectrum in numerical simulations with local non-Gaussian initial 
conditions have been presented so far, to the best of our knowledge, only in [35]. This preliminary work 
employed simulations with large values for the non-Gaussian parameter Jnl {Jnl ^ 100), pointing-out 
a peculiar dependence on Z^^, not present in the same configurations of the matter bispectrum. These 
results focused on the effects of non-Gaussianity on squeezed configurations and qualitatively tested the 
scale dependence of linear and quadratic corrections in /^vl , together with their dependence on redshift 
and halo mass threshold. 

We consider N-body simulations with local non-Gaussian initial conditions corresponding to /atl = 
±100, that is, characterized by a relatively small departure from Gaussianity, as suggested by CMB 
observations [1]. We improve over the preliminary results of [35] by presenting systematic mcasurcirncnts 
of all triangular configurations shapes from large down to mildly nonlinear scales, along the lines of our 
previous work on non-Gaussian effects on the matter bispectrum [61]. We consider a low and high mass 
halo sample and a unique output redshift at z = 0.5 (which is of the order of the median redshift of 
forthcoming redshift surveys). This allows us, in the first place, to provide an estimate of the signal 
generated by local PNG and compare it to that induced in the halo power spectrum. For all triangular 
configurations, we compare our measurements to the model of [60], where some of the scale-independent 
bias factors are fitted to the data rather than derived using the peak-background split prescriptions. A 
comparison of best-fit values for the bias parameters with the peak-background split predictions is also 
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presented. 

This paper is organized as follows. In Section II, we review the nmltivariatc bias expansion proposed 
by [25] and, in Section III, we spell out the tree-level expression for the halo bispectrum derived by 
[60] together with the analogous expression for the matter-matter-halo cross-bispectrum. In Section 
IV, we describe the numerical simulations employed in our work and, in Section V, we estimate the 
cumulative signal-to-noise for the non-Gaussian effects for both the power spectrum and bispectrum. 
In Section VI, we perform a detailed comparison between the model and the measurements of the halo 
and matter-matter-halo bispectrum, and discuss the values of the best-fit bias parameters in light of the 
peak-background split expectations. We also perform a simple Fisher matrix analysis based on the halo 
bispectrum model, in an attempt to provide a first estimate of the ability of the halo bispectrum to 
constrain PNG. We present our conclusions in Section IX. 



II. LARGE-SCALE BIAS, fNL, AND THE PEAK-BACKGROUND SPLIT 



In this section, we derive an expression for the Eulerian halo overdensity 5^ as function of the nonlinear 
and non-Gaussian matter overdensity 5 and oi the initial, i.e. linear, curvature perturbations (f) including 
corrections at second order in 5 and proportional to the product (j) 5. In doing so we follow the approach 
of [25] and [60] and assume a multivariate bias expansion of the halo overdensity. The expression for 5h 
is obtained by applying the peak-background argument along the lines of [17] and [25]. Note that we will 
not take into account the additional non-Gaussian bias corrections computed by [33, 34] and [32] since 
they are negligible for the local model with constant fNL- 

Throughout this paper we will assume non-Gaussian initial conditions of the local kind and work 
consistently at linear order in the non-Gaussian f^L parameter. This is justified in light of the limits 
provided by CMB [8] and LSS [17, 38, 39, 62] observations of a local /atl: typically |/jvi,| < 100 at 
99% CL. As we will see in Section VI, this is a good approximation even for the comparison with our 
simulations, which assume /atl = ±100. 

We shall work with a local expansion of the form [25, 60] 

5h ^bio5 + boi <Po + ^ b2o + bn (/)o (5 ©(/^l) (1) 

where S still stands for the nonlinear and non- Gaussian matter density contrast, while represents 
the linear and Gaussian curvature perturbations. This expression will be used to compute the leading, 
"tree-level" contribution to the halo bispectrum whereas, for matter and matter-curvature correlators, 
we will consider additional perturbative corrections. 



A. Lagrangian bicis 

For local quadratic non-Gaussianity, the Bardeen's curvature perturbation in the matter dominated 
era is given by 

$(X) = 0(X) + /ArL0(x)2 , (2) 

where is a Gaussian field. In the peak-background split framework, we can separate the perturbations 
into their long- wavelength and short- wavelength piece, and c^i , and thus obtain 

$(x) = </>o(x) + fNLM^f + [1 + 2/jvL<^o(x)] <^i(x) + fNLM^f + const. (3) 

In this expansion, the most relevant term is {1+2 f^hMi'i since, in a region where (fio takes some constant 
value, it can be interpreted as a local, scale-dependent rescaling of the amplitude of short-wavelength 
fluctuations. As shown in [17], the second and fourth terms can be ignored as far as one is interested in 
the fc-dependent bias correction. At linear order, the Fourier modes of the density perturbations 5 are 
related to those of the curvature perturbations $ through the linearized Poisson equation 

5{k,z) = M{k,z)^{k), (4) 

where 

(5) 
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with T{k) representing the matter transfer function and D{z) the hnear growth factor, while flm and Hq 
are the present time relative matter density and the Hubble parameter, respectively. 

Taking 0o constant and convolving the left- and right-hand side of Eq.(3) with M, we can see that the 
effect of non-Gaussian initial conditions can be interpreted as a local, /jvL-dependent modulation of the 
r.m.s. amplitude ui of the short-wavelength density perturbations (for a more rigorous derivation, see 
[28, 33]) 

cTi ^ai(l + 2/jvL0o), (6) 



Following [25, 60], we define the Lagrangian bias parameters bfj{Mi,6i) from the expansion of the La- 
grangian halo density field Sj^ in the large scale, linear matter density and curvature perturbations 

= b^odo + b^i(t>o + \b^oSo + b^iSo(t>o + .-- (7) 

Assuming universality, the bias parameters can be derived from the shape of the unconditional mass 
function /(i/) alone, provided that we substitute the variable v = Sc/(ti with 

^1 - ^0 

/I , Of — jT^ — ' 
(1 + 2fNL9o)cri 



so that, in particular 



bio 



/(i/io) dSo 



(9) 
(10) 

<5o.</'o=0 



Notice that the halo mass function itself depends on the non-Gaussian parameter /jvl- If we factorize 
the effect of non-Gaussianity on the mass function as 

f{u) = fG{l^)RNG{y). (11) 

where fai^) is the halo mass function for Gaussian initial conditions and Rng represents the relative 
effect of PNG, we can split the scale-independent bias parameters 6fo and 620 into 

= ^-fo.G + A6fo,ArG . (12) 

bw = ^20,G + ^^20,AfG ' (13) 

where the Gaussian components 6io,G are obtained from the Eq.s (9) and (10) in terms of the Gaussian 
mass function fa, while the non-Gaussian corrections are given by 



_ 1 d'^RNGM 

^''^O.iVG - 3^2 



(14) 



+ ^Ko,G^Ko,NG- (15) 



As shown in [25], all the bij with j ^ can be written in terms of the bio, so that in particular 

6^1 = 2/wi5,6fo> (16) 
6fi = 2/^^(^eb2o-6fo) • (17) 

One should keep in mind that these relations are strictly valid for a universal mass function. For non- 
universal mass functions, the non-Gaussian bias corrections should be computed through a direct evalu- 
ation of derivatives of the halo mass function with respect to mass [32] . 
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B. Eulerian hias 



The Eulerian halo density 6h can be expressed in terms of Sf^ as [63] 

Sh{S,M=S+{l + S)SUSo,^o), (18) 

where we notice now the additional dependence of Sf^ (and 5h) on the linear and Gaussian curvature 
perturbation 0o- To obtain an expression of Sh as a function of the nonlinear matter overdensity S, we 
need to express the linear density perturbations So as a function of 5. This is usually done assuming the 
expansion — J2^i ^i'^^ = ^i*^ + '^2<5^ + ••• derived in the spherical collapse approximation. In this case, 
the series coefficients correspond to spherical averages of the kernel of the perturbative expansion, with, 
in particular, a\ = 1 and 02 = —17/21^. 

On inserting Eq. (1) into Eq. (18), we can express the Eulerian bias parameters in terms of the 
Lagrangian ones [25, 60]^ 

610 = l + 6fo> (19) 
601 = 601: (20) 
h20 = 6^0 + 2 (ai + a2) > (21) 

611 = ai6n + &oi- (22) 

As in the Lagrangian case, Eq. (17) and (16), the two extra parameters induced by PNG can be written 
as a function of 610 and 620 

&01 = 2W,(6io-l), (23) 



611 = 2/, 



NL 



Scb20+ ( ^^c-l] (bio-l) 



(24) 



where we replaced the coefficients ai and a2 of the spherical collapse expansion with their numerical 
values. We note that these additional terms vanish for an unbiased population (610 = 1, 620 = 0). 

In our calculations, we will assume the above expressions for &01 a-nd bn to depend only on the Gaussian 
component of &10 and 620- In fact, in order to explicitly account for all the corrections linear in /nl, we 
split 610 and 620 into their Gaussian and non-Gaussian pieces: 

ho = ho,G + A6io,jvG , (25) 

b20 = &20,G + A620,JVG , (26) 

with 

A6io,jvG = A6fo,;vG. (27) 

A620,iVG = A62o_;vG + ^Abfo.jvG • (28) 
Taking into account Eq.s (14), (15) and (8) we have 

A?, - ^ ^ dRNG .„QX 

A620,ArG = TT^ {^^0,G- W7 ] ^ 7. — • (30) 

Si Rng ou^ 5c \ 21 J Rng ov 

where the expression for A6io,jvG has been derived in [21]. 



^ One can in principle consider the more general case where the relation between i5o and <5 is given, in Fourier space, by 
the perturbative expansion 

SqM = Sk- J d^qi d^giSoO^ - qi2) -F2(qi, q2) <5qi i5q2 + . . . . 

This can have relevant consequences on the correction to the quadratic bias due to PNG, when the Gaussian, linear 
perturbation So in term feoi^o k = f>()iSo.k/^I{k) is replaced by its the nonlinear and non-Gaussian counterpart S as 
suggested by the results of [32]. We will not further discuss this rather thorny issue, leaving it for future work. 
^ Our expression for bn is different from the same expression in [25], where bn = aibfi/2 + feg^. 
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III. THE MODEL 



As already mentioned, we will not consider any loop-correction to the halo power spectrum and bispec- 
triim induced by the nonlinear local bias expansion. Here, we will focus on the leading-order, tree-level 
expressions for the matter-matter-halo cross-bispectrum and for the halo bispectrum. They will be 
compared to the simulation results in section VI. For comparison, we will also consider the matter- 
halo cross-power spectrum and the halo power spectrum, already studied in numerical simulations in 
[9, 21, 23, 24, 29, 32]. 

The starting point is Eq. (1) which, in Fourier space, up to quadratic corrections in the density and 
curvature fields reads 

4(k) ~ 6io<5(k) + 6oi<^o(k) + ^62o ^ d^fcid^fca <5D(k - ki2) <5(ki) .^(ka) 

+ &11 J Skid^k^ (5i3(k - ki2) 0o(ki) (5(k2) + . (31) 

Here and henceforth we denote kjj...j^ = kj^ H -|- kj^. To highlight the "scale-dependent" corrections 

to halo bias factors, it is convenient to rewrite this expansion in terms of the fully nonlinear and non- 
Gaussian matter density 8 and of the Gaussian component of the imear matter density Sq. In the following 
we will assume that 5q, corresponds to the Gaussian component. We then have 

5h{\<i) ^ 6io<5(k)+coi(fc)(5o(k) + ^b2o j (i3fcid3fc2 ^^(k - ki2) <5(ki) 5(k2) 

+ ^ j d^hd^k2 ^D(k - ki2) [Cii(fci) ^o(kl) ^(k2) + (fci O k2)] , (32) 

where we have defined 

<.,W=j^, (33) 

cn(t) (34) 

coi(fc) clearly corresponds to the usual non-Gaussian scale-dependent correction to linear bias, whereas 
cii(fc) is the analogous contribution to the quadratic bias. 



A. The matter-halo and halo power spectrum 



To calculate the two-point correlators, we will neglect the contribution of loop corrections induced by 
nonlinear bias and, therefore, limit ourselves to the two linear terms of Eq. (32). 

The leading order contribution to the halo-matter cross-power spectrum is then simply given by 

P5h{k) = bwP5{k) + coi{k)Pssoik), (35) 

while for the halo power spectrum we have 

Phik) = bio Psik) + 2 6io coi(fc) PsSoik) • (36) 

Both expressions simply correspond to the usual halo bias correction of [9] plus the scale-independent 
correction introduced by [21]. Here, Pg is the nonlinear matter power spectrum whereas PgSo is defined 
as ((^kii^o.G.kg) = <5z5(ki2)f'<5<5o (A:i). We will compute these quantities in standard perturbation theory 
including loop corrections corresponding to 0{6q) in the linear density field. We refer the reader to 
Appendix A for details of the PT expressions. 

At nonlinear level, both Pg and Pssa receive correction due to the non-Gaussian initial conditions 
and, therefore, depend on /nl- As we will do later for the three-point matter correlators, we separate 
the Gaussian from the non-Gaussian contribution and write Ps = Ps,g + ^Ps,ng and Psso = PsSo,G + 
APsso,NG- The non-Gaussian corrections AP^^jvg and APsso,ng will be evaluated at linear order in Jnl- 
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Similarly, we can distinguish the Gaussian and non-Gaussian components in the matter-halo and halo 
power spectrum. For Gaussian initial conditions, the matter-halo cross-power spectrum reduces to 

P5h.G{k)=bio,GPs,Gik), (37) 
while the non-Gaussian correction is given by 

APsh,NGik) = ho,G AP6,NG{k) + Aho,NG + Coi(fc) PgSo^G + • (38) 

For the halo power spectrum we have instead 

Ph,G{k) = bloPs{k), (39) 
with the non-Gaussian correction given by 

APh,NG{k) = ?^?o,G APs,NG{k) + 2 feio,G Abw^NcPsMk) + 2 ho,G coi(fc) Pss„,G + OUIl) ■ (40) 

This distinction between Gaussian and non-Gaussian contribution will be crucial in the analysis of Sec- 
tion VI as we will compare them separately with the simulations results. 

B. The matter-matter-halo and halo bispectrum 

We define the cross matter-matter-halo bispectrum Bssh as 

((5(ki),5(k2)<5,,(k3)) EE5z5(ki23)S55^(fcl,fc2;fc3), (41) 

where we assume that the third wavenumber kz to represent the halo overdensity variable. At leading 
order in the local bias expansion, the expression for the matter-matter-halo cross-bispectrum is given by 

B55h{ki,k2;ks) = bioBs{ki,k2, k^) + 001(^3) BssSa{ki,k2\ fcs) + Ps{ki)Ps{k2) 

+ ^ [cii(/ci) Ps5o{ki)Ps{k2) + (fci ^ k2)] + 0{.f%L) , (42) 

where, in addition to the matter power spectra Ps and Pssa already introduced, Bs represents the non- 
linear matter bispectrum while Bsss„, defined as ((S(ki)(5(k2)(5o(k3)) = <5z5(ki23) -85550(^1, fc2; fc3) is the 
cross-bispectrum between two nonlinear matter density fields 5 and one linear (and Gaussian) mass over- 
density 5o. Again the last argument, ^3, refers to the 60 variable. For Gaussian initial conditions Bgsh 
reduces to 

Bssh,G{ki,k2;kz) = bio^c Bs,G{ki,k2,k3) + 620,0 ^'.5,0(^1)^5,0(^2) , (43) 
while the non-Gaussian correction explicitly is 

ABssh,NG{ki,k2; fcs) = 610,0 ABs^Noiki, fc2, ^3) + A6io,wo Bs.G{ki,k2, kz) + 001(^3) Bssso,G{ki,k2; kz) 

+b2o,G [APs^NG{ki)Ps,G{k2) + {ki ^ ka)] 

+Ab2(i,NG [Ps.,G{kl)Ps,G{k2) + {kl ^ ^2)] 

+ ^ [cii(A:i) PsSo,Giki)Ps,Gik2) + (fci ^ fca)] + 0(/^i) . (44) 
For the tree-level halo bispectrum, we obtain 

Bh{ki,k2,kz) = blQBi{ki,k2,kz) + b\^ [cai{k?,)B ssSo{ki, k2\ kz) +2^exm.] 
+6?o 620 [^^5(^1)^5(^2) + 2 perm.] 

+ ^ 610 620 [coi{ki)Psso{ki)P5{k2) + CQi{k2)PdSo{k2)Ps{ki) + 2 perm.] 
+ 1 bio [cii(fci) P56o{ki) P5{k2) + cii(fc2) Ps5o{k2) Psiki) + 2 perm.] + Oif^^) .(45) 
For Gaussian initial conditions, the halo bispectrum is 

Bh,G{ki,k2, kz) = 6?o,0 ^5,0(^1,^2, ^3) + bjo^c ^20,0 [P6,G{ki)Ps,G{k2) + 2 perm.] + 0(/^i) . (46) 
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while the non-Gaussian correction is given by 

ABh,NG{h,k2, ka) = b^ Q ABs,NG{ki,k2, fcs) + 3 Q A6io,jvg -85,0(^1, ^2, ^3) 
+blo [coi{ka)Bsss„{ki,k2; fca) + 2 perm.] 

-|-6lO,G (2 A6io,JVG ^20,G + &10,G A620,AfG) [f'5,G(fcl)-P<5,G(fe) + 2 pcim.] 

+&?o,g''2o,g [APs,NG{ki)P5,Gik2) + APs,NG{k2)Ps,G{ki) + 2 perm.] 

+ ^ &io,G&20,G [coi{ki)P56o,G{ki)P5,G{k2) + coi{k2)P66o,G{k2)P5,G{ki) + 2 perm.] 

+ 1 blo,G [cii{ki) PsSoMf^i) Ps,G{k2) + cii(A;2) P5s„G{k2) Ps,G{ki) + 2 perm.] 
+0{fl^). (47) 

We refer again the reader to Appendix A for the details of the evaluation of the matter correlators, 
including the bispectra and Bsss^, in perturbation theory. 



IV. SIMULATIONS AND HALO SAMPLES 



To measure the effect of non-Gaussian initial conditions of the local quadratic type on the halo bis- 

pcctrum, wc use a series of eight twelve sets of three N-body simulations of the ACDM cosmology, each 
of which has Jnl = 0,±100 [21]. The same Gaussian random seed field </> is employed in each triplet of 
runs so as to minimize the sampling variance. These simulations evolve 1024^ dark matter particle in a 
cubical box of size 1600 Mpc. The force resolution is 0.04 times the mean inter-particle distance. The 
(dimensionless) power spectrum of the Gaussian part </'(x) of the Bardeen potential is the usual power-law 
A^{k) = k^ P,f,{k) / = A^{k/ko)'"^~^. The spectral index is = 0.96 and the normalization of the 
Gaussian curvature perturbations is = 7.96 x 10""'^'' at the pivot point ko = 0.02Mpc~-'^, close to the 
best-fitting values inferred from CMB measurements [64]. We will consider a single output at redshift 
z = 0.509. 

Friends-of-friends (FoF) halos were extracted using a linking length of 0.2 times the mean inter- 
particle distance. We will present results for a low-mass bin defined by 8.8 x 10^^ Mq < M < 
1.6 X 10^^ M0, and a high-mass bin given by M > 1.6 x 10^^ h^^ M©. The mass thresholds has been 
chosen in order to have the same halo number density of about ri/j = 1.8 x 10"'' Mpc~^ for both mass 
bins. As we will see, the linear bias for the low mass sample is about 1.6 while for the high mass sample 
we have about 2.3. In addition the choice provides two halo populations characterized by quadratic bias 
parameters of different sign: negative for low mass halos and positive for high mass halos. This has a 
direct implication for the sign of non-Gaussian corrections due to nonlinear bias. We finally remark that 
the shot noise contribution to the power spectrum is comparable to the halo power spectrum itself at 
about 0.14 /iMpc^^ and 0.2/iMpc^^ for the low and high mass bin respectively. In all our results we 
present the halo power spectrum Ph and the halo bispectrum Bh corrected for shot-noise, while we ignore 
such correction for cross-correlations between matter and halo density fields. 



V. SIGNAL-TO-NOISE 



To assess the effect of non-Gaussianity on the halo and matter correlation functions, it is useful to 
compare the cumulative signal-to-noise for various statistics. In the case of the the power spectrum, this 
quantity is defined as 

fSVj^- [Pm,NG{k) - PrnMk)? 

\NJp Var[P™,G(fc)] ' ^ ^ 

where Pm.NG a-nd Pm.G represent the matter power spectrum for non-Gaussian and Gaussian initial 
conditions, respectively, while Ya,T[Pm,G] is the variance of the mass power spectrum in the Gaussian 
case. Given a surveyed volume ^ oc L^, the sum runs over all wavenumbers k in steps of the fundamental 
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FIG. 1: Left panel: Cumulative signal-to- noise for the effect of non-Gaussian initial conditions on the matter 
power spectrum and bispectrum as defined in Eq. (48) and Eq. (49) as a function of the largest wavenumber 
included {kmax)- Continuous and dashed curves indicate the signal-to-noise for the matter bispectrum and power 
spectrum, respectively, when the correlators and their variance are determined from the simulations. The dotted 
line correspond to the tree-level prediction for the matter bispectrum and its variance in Eulerian PT. Right 
panel: Cumulative signal-to-noise for the matter bispectrum with the sum in Eq. (49) restricted to triangles with 
fea < O.OS/iMpc"'^ (dotted curve), ks < 0.05 /iMpc^'^ (dot-dashed curve) and fes < 0.1/iMpc~^ (dashed curve) 
compared to the case where all triangles are included (continuous curve). 



mode kf = 2'k/L. Similarly, the signal-to-noise for the bispectrum is 



^^b",,<,^,3=,^ Var[B„,c(fci,fe,fc3)] • ^ ^ 

Here the sum is over all triangular configurations, i.e. all the triplets ki, k2 and ks forming a closed 
triangle. For simplicity, we only include in both cases the variance of the correlators while, in principle, the 
complete covariance between fc-bins or bispectrum triangles should be considered. Since this calculation 
is for illustrative purposes only, we just note that the effects of covariance can significantly reduce the 
signal-to-noise for both the power spectrum and bispectrum, but they marginally aflfect the comparison 
between the two correlators at large and mildly non-linear scales (see, for instance, [53, 54] for a realistic 
estimate of the effects of covariance on the galaxy power spectrum and bispectriim). 

In the left panel of Fig. 1 we show the cumulative signal-to-noise for the matter power spectrum (dashed 
curve) and bispectrum (continuous curve) as measured from the simulations. The dotted curve represents 
the predicted signal-to-noise for the matter bispectrum, assuming the correction induced by non-Gaussian 
initial conditions is given by the linearly evolved initial bispectrum Bq and the bispectrum variance is 
given by its leading Gaussian component [51, 65] 

Yar[B(k,,k2, k^)] = g^^fcffc^fcg -P(^i)-P(^2)P(fc3) , (50) 

with SB = 6, 2 or 1 for equilateral, isosceles or scalene triangles respectively. The theoretical prediction 
agrees reasonably well with the measured signal at large scales. On the other hand, the excess in the 
measured cumulative signal-to-noise at small scales, which are affected as well by a larger variance, is due 
to the impact of non-Gaussian initial conditions on the nonlinear evolution of the matter bispectrum. 
Most importantly. Fig. 1 clearly shows that the effect of PNG in the matter bispectrum is larger by roughly 
a factor of 4 than in the matter power spectrum, where PNG only enters at the nonlinear level. This is 
true on mildly nonlinear scales, but it is possible that this behavior extends to smaller scales as well [66]. 
In the right panel of Fig. 1, we show the cumulative signal-to-noise for the matter bispectrum restricted 
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FIG. 2: Cumulative signal-to-noise for the effect of non-Gaussian initial conditions as defined in Eq.s (48) and (49) 
but evaluated respectively for the halo power spectrum Ph {dashed curves) and halo bispectrum Bh {continuous 
curves), as a function of the small scale included (kmax)- Also shown are the signal-to-noise for the halo bispectrum 
restricted to triangles with kg < 0.03/iMpc~^ {dotted curves) and, for comparison, the signal-to-noise for the 
matter power spectrum and bispectrum {respectively dashed and continuous black, thtn curves). Left panels 
shows the results for low mass halos {blue curves) while right panels show the results for the high mass halos {red 
curves). Inset labels are common to both panels. 



to different classes of triangular configurations. Specifically, the sum in Eq. (49) runs over triangles 
with fca < 0.03 ft. Mpc"^ {dotted curve), < 0.05 /iMpc~^ {dot-dashed curve) and < 0.1/iMpc~^ 
{dashed curve). The dotted curve, for instance, corresponds mainly to squeezed triangles, with one 
mode (corresponding to k^) deeply in the linear regime and the two others (for which fci, ^2 > k^) 
gradually probing nonlinear scales as kmax increases. Despite the fact that squeezed configurations 
provide a significant contribution to the signal-to-noise as is expected for local non-Gaussianity, the 
signal is fairly equally distributed over all configurations. For instance, for fcmax = 0.2 ft Mpc"'^, triangles 
with ks < 0.03ftMpc^^ account for one-third of the total signal solely. 

This picture changes substantially when we consider halo correlation functions. For the halo power 
spectrum, the scale-dependent correction to the linear bias dominates the signal at large scales. In Fig. 2 
we plot the cumulative signal-to-noise defined in Eq.s (48) and (49) as estimated from measurements 
of the halo power spectrum P/j {dashed curves) and bispectrum Bh {continuous curves), respectively. 
Results are shown for the low mass {blue curves, left panel) and the high mass halo sample {red curves, 
right panel) assuming /tv^ — 100. In addition, we show the signal-to-noise for the halo bispectrum 
restricted to triangles with fcs < 0.03ftMpc~^ {dotted curves) and, for comparison, overlay the signal- 
to-noise for the matter power spectrum and bispectrum {respectively dashed and continuous black, thin 
curves) reproduced from Fig. 1. 

It is clear from all the cases considered here that the distribution of the signal as a function of scale differs 
significantly between the power spectrum and the bispectrum. In the power spectrum case, the main effect 
of non-Gaussian initial conditions is the scale-dependent correction to the linear bias. Since it is largest 
at the largest scales, the cumulative signal reaches its maximum already at fc ^ 0.05ftMpc~^, nonlinear 
corrections at smaller scales adding little signal. For the bispectrum, the cumulative signal-to-noise has 
a different dependence on the maximum wavenumber kmax mainly because it is a sum over all possible 
triangular configurations, whose number grows as fc^aa; [53]- At large scales, the signal is suppressed 
relative to that of the power spectrum since there are only a few measurable triangles characterized by a 
large variance. For larger values of kmax the number of triangles grows considerably. As a consequence, 
the effect of PNG becomes larger in the bispectrum at relatively small kmax ^ O.lftMpc"^ for the 
halo samples considered. The fact that the matter bispectrum includes the initial components linearly 
extrapolated as well as relatively larger nonlinear corrections contributes to enhance this effect, although, 
as we will see, the overall signal in the halo bispectrum is characterized as well by significant cancellations 
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of individual terms on the r.h.s. of Eq. (47). 

It is interesting to notice how the distribution of the signal among different triangles differs between 
halos and mass. As mentioned above, the dotted curves correspond to triangles with one side being less 
than 0.03/iMpc~^, i.e. mostly squeezed configurations. For halos, this subset of triangles accounts for 
a significantly larger fraction of the total signal-to-noise. As we will see later, this is related to the fact 
that the primordial component of the matter bispectrum and the effects of PNG on linear and quadratic 
bias both peak on squeezed configurations. At first sight, this is good news as one may think that most 
of the information on local non-Gaussian initial conditions can be extracted from a limited number of 
squeezed configurations. However, one should bear in mind that such triangles arc expected to be highly 
correlated as they share at least one wavenumber. In this regards, we emphasize that, even though the 
cumulative signal-to-noise shown in Fig. 1 and 2 is obtained from measurements in N-body simulation, 
it does not include the covariance among different triangular configurations. We expect it to be quite 
significant for triangles sharing one or more sides, particularly at large scales [53, 54, 67]. 

VI. RESULTS 
A. Measurements and analysis 

In this section we present measurements of the cross matter-matter-halo bispectrum Bmmh and the 
halo bispectrum in the simulations described in Section IV. To facilitate the comparison with recently 
published results, we show as well measurements of the cross matter-halo power spectrum and halo power 
spectrum. These measurements provide a first assessment of the model developed in Section III. The 
comparison between N-body measurements and predictions will involve fitting for several constant bias 
parameters such as the Gaussian component of the linear and quadratic bias 6io,g and 620,0 and their 
scale-independent, non-Gaussian corrections Abio^No and Ab2o,NG- The best fit values will be compared 
with expectations from the peak-background split theory in Section VII. 

For each of the four correlators, we will consider three specific measurements: 

1. the Gaussian component, i.e. the measurements of the correlators in the simulations with Gaus- 
sian initial conditions, that is Pmh,G, Ph,G, Bmmh,G and B/,,g; 

2. the non-Gaussian correction, i.e. the difference between the non-Gaussian (/jvl = +100) and 
the Gaussian measurements, 

ACng = CngUnl = +100) - CgUnl = 0) , 

where C stands for any correlator such as Pmh, Pm, Bmmh and Bh; this difference is first computed 
for each realization separately and then averaged over the available realizations, thereby reducing 
the scatter as each Gaussian/non-Gaussian simulations pair is obtained from the same Gaussian 

seeds; 

3. the 0{flij^) component, obtained as the mean over all realizations of the combination 

[CngUnl = +100) + CngUnl = -100) - 2Cjvg(/jvl = 0)]/2 = 0(/^i) , 

measured for each /nl = 0, ±100 triplet; this quantity is sensitive to any non-Gaussian correction 
beyond linear order in the nonlinear parameter and, therefore, is a measure of the amplitude of the 
corrections quadratic in /nl, which are neglected in our model. 

In other words, we do not compare the model directly to the measurements of each correlator, but we 
analyze separately the Gaussian and non-Gaussian components. The prediction for the non-Gaussian 
correction are easily obtained from the perturbative expressions of our model, Eq.s (38), (40), (44) and 
(47). 

We are interested, in the first place, in the behavior of the halo correlators at large scales and in the 
ability of the tree-level approximation for the bias expansion to capture the main effects of non-Gaussian 
initial conditions on the matter-matter-halo and halo bispcctra. In particular wc want to verify if the 
functional form of the different components in this approximation present the correct dependence on 
scale and on the triangle shape. In other words, we want to establish if the terms contributing to Bmmh 
and Bh, Eq.s (44) and (47) respectively, can provide an accurate model for all triangular configurations 
at large scales. 
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As mentioned above, we will treat the (scale-independent) Gaussian bias factors 610, g a-nd &20,g and 
their non-Gaussian corrections Afeio.TVG and Ab2o,NG as free parameters. While these parameters control 
the amplitude of several terms in the expressions of the halo correlators, our choice docs not limit signif- 
icantly the predictivity of the model. Firstly, we assume that the terms generating the scale-dependent 
bias corrections in the halo density expansion Eq. (32), like Coi(fc) and Cii(fc), are fully determined by 
bia,G and 620, G (which can be measured from the Gaussian realizations). Secondly, the fit involves aZ^the 
triangular configurations down to a certain scale or, equivalently, a large number of degrees of freedom. 
In our model, varying those four free parameters (as we will see shortly, there are in practice less than 
four) can affect the shape dependence of the halo bispectra only to a limited extent. 

The fitting procedure of the halo power spectrum and bispectrum follows these general steps: 

1. the Gaussian linear bias param,eter biQ,a is obtained from the Gaussian halo power spectrum Ph,G'i 

2. the non-Gaussian scale-independent correction to the linear bias Abio.NG is determined from the 
non-Gaussian correction to the halo power spectrum, AP^^ng = Ph,NG — Ph,Gt assuming the best-fit 
value of 610, G obtained from Ph,G: 

3. the Gaussian quadratic bias b-zo.G is determined from the Gaussian halo bispectrum Bh,G upon 
setting 610, G to its best-fit value from step 1. 

4. the non-Ga,ussian scale-independent correction to the quadratic bias Ab2o^NG is computed from the 
non-Gaussian correction to the halo bispectrum, ABh^Nc = Bh^NG — Bh,G^ assuming the best-fit 
values of 610, Gi A610, atg a^nd 620,0 from steps 1, 2 and 3, respectively. 

Analogously, we independently fit the model of the cross matter-halo power spectrum P^h and matter- 
matter- halo bispectrum Bmmh- 

In this way, the fit to each of the measured quantities Ph,G, AP^.ng, Bh.G and AB^^ng effectively is 
a one-parameter fit. Clearly, we could have adopted other fitting procedures where, for instance, all the 
parameters are obtained from measurements of the bispectrum alone. Keeping the model as predictive 
as possible is the main motivation of our choice. A comparison between the best fit bias parameters and 
the values predicted from the peak-background split justifies a posteriori our approach. These aspects 
will be discussed in details in Section VII. 

Since we can only expect the tree-level approximation from the bias expansion to be valid at large 
scales, wc restrict the fits to wavcnumbcrs k < O.OTft.Mpc""'^ for both the power spectra and the three 
sides of the triangular bispectrum configurations. This choice notwithstanding, wc will show in the figures 
measurements of the power spectrum and bispectrum up to 0.2 /iMpc^^, together with the extrapolation 
of the theoretical model. Therefore, one should keep in mind that the predictions shown are not obtained 
from a fit to all the data points displayed in the figure. 

In addition, it is important to note that the fits to the bispectra, although limited to large scales, assume 
all measurable configurations up to the aforementioned kmax = 0.07 h Mpc~^ . The bispectra are measured 
for triangle sides which are multiples of Ak = 3kf = 0.012 ft,Mpc^^, where kf = 2tt/L ~ 0.004 /iMpc""^ is 
the fundamental frequency of the box. It follows that, for the halo bispectrum, we measure 597 triangular 
configurations characterized by /cs < ^2 < fci < 0.2/iMpc~^ and k^ > ki — k2, which ensure no double 
counting and enforce the triangle constraint. However, only 43 out of the 597 triangle configurations, 
for which ^3 < ^2 < fci < 0.07 /i Mpc"'^, are eventually considered for the fit of the bias parameters. 
For the matter-matter-halo bispectrum, the variable fcs is corresponds to the halo density wavemode. 
Therefore, its value is allowed to vary in the range fci — ^2 < fcs < fci + k2 since, in this case, values of k^, 
greater than k2 do not correspond to double counting. The total number of configurations measured for 
the matter-mattcr-halo is 1,549, among which 93 are used for the fits. While in figures 6 to 13 wc will 
present results for specific bispectrum configurations only, wc stress once more that wc arc not simply 
fitting the triangles shown at any given time but perform instead a global fit to all large scales triangles 
of any shape. 

B. Power spectra 

In Fig. 3 and 4, we show results concerning the cross matter-halo power spectrum P^h and the halo 
power spectrum Ph, respectively. Each figure displays the measurements of the Gaussian correlator Pq 

and the residuals to the model fit {first two rows), the non-Gaussian correction AP^vG and the relative 
residuals {third and fourth row), the ratio non-Gaussian to Gaussian P^g/Pg and the measurement of 
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FIG. 3: Cross matter-halo power spectrum, Psh{k). See text for explanation. 




FIG. 4: Halo power spectrum, Ph(k). See text for explanation. 
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the 0{f'j^]^) component {last row). Left panels correspond to the low mass bin 8.8 x lO^^ft-^^Mo < 
M < 1.6 X 10^^h~^MQ (data points in blue), while right panels correspond to the high mass bin 
M > 1.6 X 10^'^/i~^Mq (data points in red). The shaded gray area in the residual plots indicate a 
1% and 10% error for the Gaussian power spectrum and the non-Gaussian correction, respectively. The 
thin vertical line indicate the maximum wavenumber, kmax = 0.07/iMpc~^, used for the fits. 

Let us consider first the matter-matter-halo power spectrum, Fig. 3. In the upper panels, the Gaussian 
measurements are compared to the simple theoretical prediction Eq. (37), where the linear bias parameter 
6io,G is fitted for and the nonlinear matter power spectrum is computed up to one-loop in PT. Since the 
fit is restricted to fc < 0.07/iMpc~^, such nonlinear corrections are not affecting the determination of 
the linear bias. At the same time, they do not improve much the agreement between the model and the 
simulations at smaller scales, where further corrections due to nonlinear bias should be significant. At 
large fc, the model overestimates the data points for the low mass halos, and underestimates them for the 
large mass halos. As we will see shortly, this is consistent the Gaussian quadratic bias parameter being 
negative and positive, respectively. 

In the second and third rows, the correction induced by PNG, APsh,NGy is compared to the theoretical 
prediction Eq. (38) {continuous curves). The various curves represent the contribution generated by the 
scale-dependent bias correction coi(fc) {dashed curve), by the scale- independent bias correction A6io,jvG 
{dot-dashed curves: the absolute value is shown as this term is negative) and by the non-Gaussian correc- 
tions to the matter power spectrum AP^ jvg {dotted curves). As already remarked in [21] and [25], adding 
the scale-independent correction A6io,jvg to the scale-dependent bias coi{k) significantly improves the 
agreement with the simulations. At smaller scales, further improvements can be achieved on including 
the additional non-Gaussian corrections derived in the peak-background split approach of [32, 33]. Nev- 
ertheless, we shall ignore them here as we focus on the simplest model for the large-scale halo bispcctrum. 
The fourth row shows the non-Gaussian to Gaussian ratio Psh,NG/Psh,G with the same labeling for the 
various model components. Note that the contribution from the Abio,NG term is negative. 

The lower panels show the combination [Psh.NGtfNL = +W0) + Psh.NG{.fNL = ~W0)~2Psh,G{.fNL = 
0)]/2, corresponding to 0{f'^j^) contributions the matter-halo power spectrum. These are consistent with 
zero for low mass halos whereas, for the high mass halos, there is some evidence for a signal at large 
scales. These second-order corrections /nl are negative, of the order of 5% of the whole non-Gaussian 
correction APs^ nq. They might be due to our approximation for Coi(fc) ~ (6io — 1) ~ (&io,G — !)• The 
inclusion of the scale-independent correction on the linear bias parameter, Abio^NG in the expression for 
coi{k), Eq. (23), leads to a second order effect in Jnl which we neglect here and do not explore further. 

Fig. 4 shows the analogous quantities for the halo power spectrum. Considerations similar to those 
above apply to these figures. However, an interesting difference is the effect of nonlinear bias corrections 
to the Gaussian halo bias at small scales, that presents the same sign for both low and high halo masses. 
In addition, there is evidence for a fff^ correction in the lowest k bin for the high mass and, to a lesser 
extent, the low mass halos. Still, on intermediate scales 0.02/iMpc~^ < k < 0.05, the effect appears to 
have a negative sign, like in the matter-halo power spectrum. Overall, it is at most 1% of the Gaussian 
halo power spectrum Ph,G, and at most 10% of the total non-Gaussian correction APh,NG- Again, its 
inclusion might improve marginally the agreement of model with the measured AP^^ng- 

Our power spectrum measurements agree broadly with those of previous studies [9, 21, 23, 24, 29, 32, 
34]. Our choice to fit for the non-Gaussian correction to the linear bias A6io,ArG can account, among 
others, for the additional large-scale corrections predicted by the approach of [32, 34] on top of the 
standard contribution, coi(fc). Notice that, unlike [29] for instance, we do not add an extra fudge factor 
to the amplitude of the scale-dependent correction. Fitting the non-Gaussian corrections APsh.ng and 
APh,NG for both A6io_jvG and this additional parameter would lead to results consistent with those of 
[29], yet significantly complicate the bispectrum analysis. We postpone a more detailed analysis to future 
work. 

Overall the results of this section show that the linear bias models for Psh,NG and Ph,NG agree at the 
percent level with the simulations results up to fc ^ 0.07h/Mpc. The non-Gaussian correction, in fact, 
represent, over the range of scales probed, a 10% correction at most, described by the simple model for 
^Psh,NG and AP/j jvG with an accuracy of a few percent. 

C. Bispectra 

We now test the models for the cross matter-matter-halo and halo bispectra. As explained above, 
for the analysis of the Gaussian Bssh,G we assume the best-fit value for the Gaussian linear bias 6io,g 
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obtained, in turn, from the analysis of the matter-halo power spectrum. We therefore determine from 
Bssh,G only the quadratic bias parameter 620.0, which we then use, together with 610, g and A610, atg 
from Pgf^, in the fit of the non-Gaussian correction ABssh,NG- Again, the latter is a one-parameter fit 
in terms of Afe20,JVG- The same applies, independently, to the halo bispectrum components Bh^c and 

ABh,NG- 

The fits to bispectrum measurements are performed for all the triangular configurations with sides 
smaller than or equal to a given k^ax- Therefore, in order to provide a complete picture of the goodness 
of the fit, we show in the first place the value as a function of kmax- This will allow us, in particular, 
to justify our choice for kmax 

= 0.07 /iMpc ^ as the wavenumber below which we expect our theoretical 

description to be accurate. 

The as a function of k^ax is defined as 



where the variance is defined as the square of the error on the mean measured from the twelve available 
realizations, and where the sum runs over all triangles defined by wavenumbers below or equal to k^ax- 
We consider values of kmax ranging from AA; = 3A;/ = 0.012 /iMpc~^ up to 0.2/iMpc~^. The best-fit 
bias parameters are determined from configurations with sides less than kmax solely. For the sake of 
consistency, the values for the linear bias parameters 610, g and A6io,jvg are determined from power 
spectra measurements at wavenumbers up to k„,,aTi and similarly for the value of &20,G in the analysis of 
^Bssh,NG and AB^^ng- Notice that the number of degrees of freedom grows considerably as a function 
of kmax, as does the number of triangles included. 

In Fig. 5 we show the per degree of freedom as a function of k„iax for both Bssh and Bh {thick con- 
tinuous curves). More precisely, the first two rows show the x^ corresponding to the Gaussian component 
of the cross-bispectrum Bssh,G and to the non-Gaussian correction ABssh,NG- The same quantities for 
the halo bispectrum components B^^g and ABh^^G are shown in the third and fourth row. Panels on the 
left correspond to the low mass halos, while panels on the right to the high mass halos. The thin curves 
represent the same quantities obtained assuming the tree-level predictions for the matter bispectra and 
linear matter power spectra. 

Firstly, it should be noted that the x^ per degree of freedom for the Gaussian components Bssh,G 
and Bh G is roughly constant, betwecni values of 1 to 2 up to scales ^0.1 /?, Mpc^ . This indicates that, 
on large scales, the tree-level approximation from the bias expansion captures well the scale and shape 
dependence of the halo bispectrum measurements for both low and high mass halos, thereby justifying 
our choice of fitting the bias parameters using all triangles formed by wavenumbers equal or smaller 
than 0.07/iMpc~^. If smaller scales or larger wavenumbers are included, then this simple approximation 
breaks down, but it may be improved by adding further nonlinear corrections in the bias expansion or 
non-local terms. However, this problem appears to be more severe for the Gaussian components than for 
the non-Gaussian corrections ABssh,NG and ABh,NG- 

Secondly, a comparison of the thick and thin lines suggests that the inclusion of one-loop corrections 
to the matter correlators does not improves the fits significantly at large scales {k < 0.07/iMpc~^). The 
effect at smaller scales cannot be properly assessed as nonlinearities in the bias expansion, which clearly 
play a major role, arc not accounted for. 

In the following figures, we compare directly the measured bispectra with the model for a choice of 
specific triangular configurations. The chosen subsets of triangles are defined as follows 

• squeezed configurations: ^3 is kept constant at ks = Ak = 0.012 h Mpc~^ (the smallest measured 

wavcrmmbcr) while fci = A;2 = k vary from AA" to 0.02 /i Mpc""'^; the first data point corresponds 
therefore to a large-scale triangular configuration while as k grows we obtain increasingly squeezed 
configurations; 

• generic configurations (I): ki and fc2 are kept constant at ki = 0.071 /iMpc~^ and A;2 = 
0.082 ft, Mpc""'^ while the angle 6 between ki and k2 varies from to tt; this corresponds to A3 
taking values from ki + k2 = 0.15 ft Mpc~^ to ki — ki = Ak = 0.012 ft Mpc~^ {left to right in the 
plot, as 6 increases); in this case, the first data point corresponds to a folded triangle (A;3 = A;i -|- A;2) 




(51) 
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Halo bispectrum, Bh 




FIG. 5; Reduced for the matter-matter-halo bispectrum fits. Thin lines correspond to the same quantity 
obtained assuming a tree-level approximation for matter correlators. 



at relatively small scales while 9 ^ tt corresponds to a squeezed limit'^; 

• generic configurations (II): fci and fc2 are kept constant at fci — 0.047 ft- Mpc~^ and ^2 — 
0.071 ft Mpc^^ while the angle between ki and k2 varies from to tt; this corresponds to ky, 



^ Nota bene: the angle 8 is defined in terms of tfie vectors ki, k2 and ka satisfying the property ki -|- k2 -I- ka = 0, not as 
an internal angle of tfic triangle with sides ki, k2 and k^. 
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taking values from /ci + ^2 = 0.12 /iMpc ^ to ki — k2 = Ak = 0.024 ft, Mpc ^] these triangles are 
almost all scalene and approximately equally distant from equilateral or squeezed configurations; 

• equilateral configurations: ki = k2 = k^ = k with k varying from AA; = 0.012 ft, Mpc~^ to 
0.02ftMpc-^ 

These configurations are chosen to given a fair assessment of the effects of non-Gaussianity on the halo 
bispectrum which, as we will see, is not limited to squeezed configurations. 

Figures 6, 7, 8 and 9 show, respectively, the four different subsets of triangles described above for the 
matter-matter-halo bispectrum Bssh{ki,k2; k^). Notice that the third variable ks corresponds to the halo 
density contrast Sh{^3)- 

As for the power spectrum plots, for each set of configurations we show the Gaussian component Bs5h,G 
and the model residuals (upper two rows), the non-Gaussian correction ABssh.NG and the corresponding 
residuals {third and fourth row), the non-Gaussian to Gaussian ratio Bssh.Na/Bssh.G {fifth row) and the 
^Unl) component determined from the average of [Bssh,NG{fNL = +100) + Bssh,NG{fNL = -100)] - 
Bs5h,G{fNL = 0) measured in each realization {last row). Left and right columns show results for the low 
and high-mass halos, respectively. 

For the Gaussian piece of the cross bispectrum, the model is given in terms of two components, the first 
and second terms on the r.h.s. of Eq. (43), shown in the plots as a dashed and dotted curve, respectively. 
These terms are simply labelled by bio and &20) since the latter are the bias parameters controlling their 
amplitude. We implicitly assign to 610 and 620 the value corresponding to Gaussian initial conditions, 
^•C- bio^G and 620,0- For the low mass halos, &20.G is negative as is the corresponding contribution to the 
bispectrum. For this reason, we show its absolute value in the log-log plots of squeezed and equilateral 
configurations. When a given contribution is negative, we denote it in the plot legend by a negative 
sign in front of the related symbol. In the plots of the models residuals for the Gaussian term Bssh,G, 
the shaded area indicate a 10% deviation. The thin, gray vertical line at fc = 0.07 shows the subset of 
data points being part of the larger subset used for the model fit. Notice that, in the case of the generic 
configurations (I) for which ^2 = 0.082 ft Mpc~^, none of the triangles shown is used for the fit. 

For nearly all triangles, with the exception of the equilateral ones, the model and the data generally 
agree within 10% at large scales. For squeezed configurations in particular, such an agreement persists 
well beyond kmax = 

0.07ftMpc"\ 

For the non-Gaussian correction ABssh,NG we have several distinct components. The first three terms 
on the r.h.s. of Eq. (44) are shown in the plots respectively as dashed, dotted and dot-dashed black 
curves, and labelled as 610, A610 and Coi(fc) since they loosely correspond to the /wear halo bias and its 
corrections. The remaining three terms on the r.h.s. of Eq. (44), which are corrections related to the 
quadratic bias, are shown as dashed, dotted and dot-dashed magenta (or light gray) curves and labelled 
as 620, A620 and cii{k). The continuous black curve shows the sum of all the terms. In the residual plots, 
the shaded area for ABssh^^G indicate a less than 50% deviation between the model and the data. The 
same notation is assumed for the various curves in the ratio plots on the fifth row. 

Even for the correction ABssh,NG in the case of the squeezed configurations shown in Fig. 6, the 
agreement between the model and simulations extends to relatively small scales. More interestingly, 
the non-Gaussian signal in the squeezed limit results from a comparable effects of PNG on the matter 
bispectrum and on the linear bias, coi(fc). For these specific configurations, the effect on the quadratic 
bias is very small and even vanishes in the squeezed limit, because the small, constant side of the triangle 
is in this case k^, the wavenumber corresponding to the halo ovcrdcnsity, while the scale-dependent 
corrections cii(A:i) and cii(A;2) are suppressed for large values of ki and k^- 

Figures 7 and 8 confirm that non-Gaussian corrections, now smaller for generic configurations, are the 
results of two distinct contributions which are roughly of the same order. Since for these configurations 
ki and k^ are fixed, the non-Gaussian correction to the quadratic bias is constant and, in this case, 
negative. As a consequence, only the primordial component to the matter bispectrum and the correction 
to the linear bias contribute to the non-Gaussian effect in the squeezed limit, which is attained for the 
triangles with — > tt. This is particularly obvious in figure 7 where ki and fc2 take similar values. The 
constant correction to 62 is most evident for nearly equilateral configurations, especially in figure 8 where 
it amounts to a reduction of the overall non-Gaussian correction to the cross-bispectrum. What is more 
interesting is the fact that, since we are only fitting for the value of Afe2o,AfGi the model correctly predicts 
the shape-dependence of ABssh,NG- This is not trivial, since such dependency is given by the peculiar 
combination of the two terms mentioned above. 

No significant correction beyond linear order in f^L to the cross bispectrum is detected for generic 
triangles and for all halos, with the exception of high mass halos in the squeezed configurations. However, 
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Squeezed configurations, Bssh{Ak, k, k) 
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FIG. 6: Squeezed configurations of t tie cross matter-matter-halo bispectrum, B5Sh{Ak,k,k). See text for expla- 
nation. 
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FIG. 7: Generic configurations (I) of the cross-bispectrum, Bssh{ki,k2,6), with fci — 0.07 /iMpc ^, = 
0.08/iMpc~^. See text for explanation. 
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FIG. 8: Generic configurations (II) of the cross-bispectrum, Bssh{ki,k2,d), with fci — 0.05/iMpc ^, = 
0.07/iMpc~^. See text for explanation. 
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FIG. 9: Equilateral configurations of the cross matter-matter-halo bispectrum, Bssh{k,k,k) as a function of k. 
See text for explanation. 
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it is small and does not affect the analysis of the overall non-Gaussian correction in terms of a model 
linear in Jnl- 

Figures 10, 11, 12 and 13 show the same subsets of triangular configurations for the halo bispectrum 
Bh. Most of the comments made for the cross-bispectrum Bgsh can be repeated here. For Gaussian 
initial conditions, the simple tree-level, local bias prescription provides a model accurate at the 10% 
level for generic triangles, the largest deviations occurring for equilateral configurations. In the halo 
bias case, however, the measurements are considerably noisier than those of the cross matter-matter-halo 
bispectrum. 

The notation adopted to identify each term in the expressions for the Gaussian and non-Gaussian 
components on the plots is the same as in the previous figures, although additional dependencies on the 
Gaussian linear bias parameter b^o are now present. For Gaussian initial conditions for instance, the 
''lo G^20,G Ps.ciki) Ps.aiM) + eye. contribution is simply denoted by 620 in the legend of the plot. There 
is one non-trivial additional contribution in the non-Gaussian correction ABh^Nd given by the sixth term 
on the r.h.s. of Eq. (47). It is shown as a continuous magenta (or light gray) curve and labelled 620 coi{k) 
in the legends, as it depends both on the scale-dependent linear bias correction and the quadratic bias. 

We remark in the first place that, in our model, the overall non-Gaussian correction to the halo bias 
^Bfi,NG is the sum of several contribution of nearly equal importance. This is particularly evident in 
the case of squeezed triangles and for the low mass halos (Fig. 10, left column), where all quadratic bias 
terms are negative and the overall signal thus is the result of large cancellations (with the caveat that 
for these configurations and mass bin in particular the model shows a relatively larger discrepancy w.r.t. 
simulations than for other triangles). For the same triangular configurations and for the high mass halos, 
only the A620 term is negative and the largest corrections are due to the matter bispectrum and the 
scale-dependent linear bias correction, the quadratic term proportional to cii(fc) is also quite relevant, 
particularly at the largest scales probed. We notice that the model under-predicts the measurements, in 
the large mass bin, by about 20%, as is also evident from the non-Gaussian to Gaussian ratio plots. At 
the same time, however, it over-predicts by roughly the same amount the measured halo bispectrum for 
triangles of sides {A/c, k — Ak, k} (not shown), very close to the squeezed configurations considered here. 
These differences might be due to binning effects to be explored in future works. A better agreement 
between the model and simulations is found for nearly all generic configurations as can be seen in Fig. 11 
and 12. The plot of the non-Gaussian correction AB^^ng for the large mass halos {right column, third 
row of Fig. 12) is particularly illustrative of the situation, where almost all contributions have a similar 
absolute values in the — >■ tt squeezed limit and several have a negative sign. These results suggest that 
it is quite difficult to describe the effects of PNG in terms of few basic corrections, or in terms of the 
scale-dependent corrections coi(fc) and cii{k) alone. 

Furthermore, the measurements of the halo bispectrum at low masses present a significant noise, par- 
ticularly evident for equilateral configurations. We checked that such scatter is not due to any individual 
faulty realization, but seems to be proper to such halo population. 

Finally, note that 0{fff^) effects are present in the halo bispectrum squeezed configurations. In 
the high mass bin, a positive 0{f'^j^) correction is measured for triangles given by {Ak,k,k} with 
k ~ 0.025 /iMpc~^, accounting for about 20% the overall non-Gaussian effect AB^^ng- Such a correction 
is, however, not present in the small mass bin for the same triangles. 

VII. BIAS PARAMETERS 

So far we have not discussed the best-fit values obtained for the Gaussian linear and quadratic bias 
parameters 610, g and 620,0 and their non-Gaussian, scale-independent corrections Abio^NG and Ab2o,NG- 
As explained in Section VIA, the linear bias parameter 610, g is determined from the power spectrum 
measurements in simulations with Gaussian initial conditions while its correction A6io,jvg is obtained 
from the extra contribution to the power spectrum induced by non-Gaussianity. The quadratic parameter 
620, G is then given by fitting the bispectrum with Gaussian initial conditions, while the best-fit value 
of Ab2o,NG is obtained from the non-Gaussian correction to the bispectrum. The whole procedure is 
applied independently to the matter-halo cross correlators and to the halo correlators (power spectrum 
and bispectrum). 

While being likely the most "predictive" procedure which does not involve a direct evaluation of the bias 
parameters, this is by no means the only possible one. We did also consider alternative determinations 
entirely based on bispectrum measurements. The outcome of these different procedures is shown in 
Fig. 14 and compared to the theoretical predictions of the peak-background split approach [63, 69-71]. 
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Squeezed configurations, Bh{^k, k, k) 
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FIG. 10: Equilateral configurations of the halo bispectrum, Bh{Sk, k, k), compared with the theoretical prediction 
assuming the best fit values for the bias parameters 6i and 62 • High mass bin (right). 
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Generic configurations (I), Bh{ki,k2 
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FIG. 11: Equilateral configurations of the halo bispectrum, Bh{k, k, k), compared with the theoretical prediction 
assuming the best fit values for the bias parameters bi and 62 • 
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Generic configurations (II), Bh{ki,k2,0), fci = 0.05 /i Mpc"\ k2 = 0.07 ft Mpc~^ 
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FIG. 12: Equilateral configurations of the halo bispectrum, Bh{k, k, k), compared with the theoretical prediction 
assuming the best fit values for the bias parameters bi and 62 • 
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FIG. 13: Equilateral configurations of the halo bispectrum, Bhik, k, k), compared with the theoretical prediction 
assuming the best fit values for the bias parameters 61 and 62 • 



28 




FIG. 14: Best-fit bias parameters versus theoretical predictions. In all panels the continuous curve shows the pre- 
dictions assuming a Sheth-Tormen mass function with the original parameters and the non-Gaussian correction to 
the mass function with the form proposed by [68] . The circles correspond to the best-fit bias parameters where the 
Gaussian linear bias 6io,g and its non-Gaussian, scale- independent correction A6io,jvg are determined from power 
spectrum measurements while only the Gaussian quadratic bias 620,0 and its non-Gaussian correction Ab2o,NG 
are determined from the bispectrum. The square data points correspond instead to the same bias parameters 
determined exclusively from bispectrum measurements. Filled symbols are derived from halo correlators, empty 
symbols from matter-halo cross-correlators. Data points are plotted at the mean mass value for the corresponding 
mass bin and are slightly displaced for clarity when needed. Vertical thin gray lines correspond to the thresholds 
defining the two mass bins. 



In the upper left panel of Fig. 14 the best-fit values of 610, g ^.Te shown as circles when they are 
obtained from the power spectrum measurements, and as squares when they are obtained from the 
bispectrum. These values are plotted at the mean mass for each of the two mass bin considered. Here, 
like in the other panels, filled symbols refer to halo correlators while empty symbols to matter-halo 
cross-correlators. Errors on the bias parameters are not shown, because an analysis including the power 
spectrum and bispectrum variance alone and neglecting covariances underestimate them significantly. 
Clearly, the best-fit values obtained from the power spectrum are about 10% larger than those obtained 
from the bispectrum. A similar discrepancy has been recently reported in [72] (who studied Fourier space 
correlators) and in [73] (who studied configuration space correlation functions). Overall, our findings are 
consistent with the results of these studies. While effects due to smoothing, particularly in relation to 
the scatter between and S in position space, and the limitations due to the tree-level model assumed 
could plausibly be invoked to explain such a discrepancy, the most likely explanation is the existence 
of nonlocal terms usually neglected in the halo bispectrum expression, present also for Gaussian initial 
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conditions [74]. Properly addressing this issue is, however, beyond the scope of this work. The figures 
also display the predicted value of the linear bias parameter obtained from the peak-background split 
according to Eq.s (9) and (19) assuming the Sheth-Tormen (ST) [70] unconditional mass function. 

The upper right panel of Fig. 14 shows instead the best-fit values of the quadratic bias parameter 620,0; 
obtained from the bispectrum measurements. We consider here as well two procedures. In the first one, 
which is used throughout this paper, the value of 610, g is determined from the power spectrum while 
the bispectrum only provides 620,0; the results are denoted by circles in the plot. In the second one, 
both 610, G and 620,0 are determined from the bispectrum; the results are denoted by squares in the plot. 
Clearly, the discrepancy discussed above in the determination of 610. g induces different values of 620.0- 
Computing again the quadratic halo bias from the mass function by means of Eq.s (10) and (21), the 
values obtained from the ST mass function are in good qualitative agreement with the best-fit values. 

In the lower left panel of Fig. 14 we show the non-Gaussian correction to the linear bias Abio.AfG- 
Circles indicate the values obtained from APsfi,NG and APh^NG (empty and filled, respectively). Squares 
indicate instead the values determined from ABssh,NG and ABh^NG [empty and filled). The value for the 
Gaussian component 610,0 is provided by measurements of the Gaussian power spectrum in the first case, 
and measurements of the bispectrum in the second case. The values obtained from the power spectrum 
are in relatively good agreement with the theoretical prediction of Eq. (29), which is obtained from the 
ratio Rng{'^) = fNG{i^)/fG{'^) of the non-Gaussian to Gaussian mass function (see Section II). For this 
quantity, we use the expression of [68] based on an Edgeworth expansion of the Press-Schechter Gaussian 
mass function, [75]. In what follows, we limit the expansion to include linear corrections in f^L only. We 
consider an expression given by Rng{qi^) where we fit for the shift parameter q comparing the prediction 
with the ratio Rng measured in our simulations, finding the best-fit value q ~ 0.91 (and very close values 
at different redshift). The non-Gaussian correction to the linear bias obtained from the bispectrum alone 
is in stark disagreement with this prediction. This, in part, motivates our choice to assume the values of 
A6io,jvo determined from power spectrum measurements. Clearly, in principle, we would obtained the 
same results using directly the theoretical predictions. 

Finally, the lower right panel of Fig. 14 shows the non-Gaussian correction to the quadratic bias 
A62o,Afo- Here as well we consider two different fits. The first one, shown by the circles in the plot and 
assumed in the previous section, consists in fitting the measurements of the non-Gaussian correction to 
the bispectrum only for A62o,jvo where the values of 610,0 and A6io,ArG are determined from the power 
spectrum and 620.0 from the bispectrum with Gaussian initial conditions. In the second case A6io.ArG and 
A62o,Aro are both determined by the non-Gaussian correction of the bispectrum, and shown by squares. 
We find again that the results in the first case are closer to the theoretical prediction, Eq. (30), with the 
exception of the result for the large mass bin from the matter-matter-halo cross-bispcctrum. 

Lastly, we note that the as a function of kmax for the comparison of the model to the bispectrum 
measurements shown in Fig. 5 does not significantly change for the different choices of parameters fitting 
discussed in this section. 



VIII. FISHER MATRIX 



The model we have tested can be used to make a preliminary forecast for the ability of galaxy bispectrum 
measurements in future redshift surveys to constrain a non-Gaussian component in the initial conditions. 
Here, we perform a Fisher matrix analysis in terms of the non-Gaussian parameter /^vl and of the 
Gaussian bias parameters 610. g and 620,0- In this section we denote the Gaussian parts of the galaxy bias 
parameters as 61 and 62, and their non-Gaussian scale- independent corrections as A61 and A62. 

This analysis is purely for illustrative purposes. It is intended to compare constraints from power 
spectrum and bispectrum measurements, but is not meant to provide specific forecasts for any future 
survey. A more detailed and realistic study will be considered elsewhere. The Fisher matrix for the 
galaxy bispectrum is thus simply defined as 



p ^ X ^ dBg{ki,k2, fca) dBg{ki,k2,k3) 1 

^ dPa dV& ABl{kr,k2M' ^ ' 



where the indices a and /3 run over the three parameters Jnl, 61, 62 while, for simplicity, we fix the 
cosmology to be that of the simulations (see Sec. IV). The Fisher Matrix for the galaxy power spectrum 

is defined in an analogous way. The fiducial values of the galaxy bias factors are obtained by computing 
integrals of the halo bias functions times the mass functions above a certain mass threshold chosen 
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FIG. 15: One-(j uncertainty on the Jnl parameter, marginalized over the Gaussian bias parameters, obtained 
from the Fisher matrix analysis of the power spectrum {dotted curve), bispectrum (dashed curve) and combined 
{continuous curve) for an ideal survey of 10 h~'^ Gpc^ at redshift z = 1, assuming a fiducial values for the non- 
Gaussian and bias parameters given by /atz, = 10, &io,g = 2 and 620,0 = 0.8. The left panel shows A/jvz, as a 
function of the maximum wavenumber included, kmax, while the right panel assumes kmax = 0.1/iMpc~^ and 
limits the smallest wavenumber included by kmin ■ 



to provide a given galaxy number density of fig. This is equivalent to an Halo Occupation Distribution 
assigning one galaxy per halo above the threshold. In addition, we assume that the scale-independent bias 
corrections A6, are fimctions of the Gaussian parameters themselves, i.e. Abi{bi), obtained by varying the 
galaxy number density. This assumption is partially justified by the strong correlation expected between 
these parameters as we vary the characteristics of the galaxy population. 

Notice that we account only for the variance of the galaxy bispectrum, given by 

AS2(fci, fc2, fcs) = ^^i^^^^Ptot{ki)Ptot{k2)Ptot{k^) , (53) 

with Ptot{k) — [Pg{k) + l/[(27r)'^ng] is the total galaxy power spectrum, including shot noise. The 
expression for the galaxy power spectrum is given by Eq.s (39) and (40) while for the galaxy bispectrum 
by Eq.s (46) and (47) where the bias parameters arc now to be interpreted as galaxy bias. For simplicity, 
we evaluate all the matter correlators at linear and tree-level for the power spectrum and bispectrum, 
respectively. In addition, the computation of the galaxy bispectrum variance is linearized with respect 
to /nl- 

The results of the Fisher analysis can be read off in Fig. 15, which shows the one-cr error on the 
non-Gaussian parameter f^i^ obtained upon marginalizing over the two bias parameter. Errors are 
plotted as a function of the maximum wavcmimbcr kmax for an ideal survey of 10 ft.^"^ Gpc^ at redshift 
z = 1. We consider a galaxy population characterized by number density fig = 10~^ Mpc~^ and by the 
Gaussian bias parameters 61 = 2 and 62 = 0.8. The dotted curve represents the error obtained from an 
analysis of the galaxy power spectrum only, the dashed curve corresponds to the galaxy bispectrum only, 
and the continuous curve is the constraint from a combined analysis of the galaxy power spectrum and 
bispectrum. As expected from signal-to-noise considerations for the effect of primordial non-Gaussianity 
on halo correlators (see Section V), the determination of /jvl from power spectrum measurements does 
not improve significantly beyond the largest scales accessible. However, due to the 

increase in the number of triangles included in the analysis, the bispectrum provides a comparable error 
A/jvL for relative small values of kmax, even before the mildly nonlinear regime. More interestingly, the 
combined power spectrum and bispectrum analysis improves that based on the power spectrum alone 
already at very large scales, even after the marginalization over the bias parameters. Wc can compare 
these results with those of [52], where the sole effect of non-Gaussian initial conditions is on the matter 
bispectrum and the non-Gaussian galaxy bias is not taken into account. We find a difference between 
the two analysis of a factor slightly larger than three, essentially due to the effect of PNG on halo bias. 

As we have seen, most of the signal in power spectrum measurements resides in the smallest wavenumber 
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available. In the analysis of observational data, it is thus crucial to avoid any systematic error that may 
arise from an improper determination of the galaxy selection function. Regarding the bispectrum, since 
the signal is distributed over a large number of triangular configurations, we can ask ourselves how 
the error on /jvz, depends on the largest scale included in the analysis, defined by the value of kmin- 
On the right panel of Fig. 15 we present the same quantities as the left panel, now as a function of 
kmin and a fixed k^ax — 0.1/iMpc~^. Clearly, combining power spectrum and bispectrum can provide 
errors comparable to those obtained from an analysis of the power spectrum alone even if the lowest 
wavemodes are excluded. This indicates that the galaxy bispectrum can provide, at the very least, a 
crucial cross-check to any power spectrum results. 

We emphasize that the Fisher matrix results presented in this section do not account for several 
important issues affecting analyses of galax;y survey data. In particular, we are neglecting the effects 
of the survey selection function and of the covariance properties of power spectrum, bispectrum and 
the cross-covariance between the two correlators. These effects are indeed responsible for a significant 
degradation of the available signal (see for instance the analysis of [54] in the context of cosmological 
parameters). Their inclusion is essential to provide realistic forecasts for any upcoming mission. This 
will be the subject of future work. 

IX. CONCLUSIONS 

In this work, we have presented the first detailed analysis of the effects of non-Gaussian initial conditions 
of the local kind on the bispectrum of halos extracted from numerical simulations. We have measured all 
triangular configurations at large scales for two different halo populations. 

We have shown that the cumulative signal-to-noise in the bispectrum exceeds the signal-to-noise in the 
power spectrum when all triangles down to mildly nonlinear scales are taken into account. The effects of 
local non-Gaussianity on the halo power spectrum are mainly due to the scale-dependent corrections to 
the linear halo bias, concentrating the signal in the smallest wavenumbers (i.e. largest scales) accessible in 
the simulations. On the other hand, non-Gaussian initial conditions of the local kind have more complex 
effects on the halo bispectrum. This leads to a non-Gaussian signal distributed across a wide variety of 
triangular configurations (with very different scales and shapes). In fact, the halo bispectrum includes, in 
the first place, the linearly evolved primordial contribution to the matter bispectrum, which is significant 
at large scales for generic models of non-Gaussianity. Furthermore, for models of PNG characterized by 
a large primordial bispectrum in the squeezed limit, scale-dependent corrections to the bias are present 
both at linear and quadratic level. As we have seen, for local non-Gaussianity these are as important as 
the primordial contribution at large scales. In general however, even scale- independent corrections are 
relevant and must be properly modeled in order to reproduce the simulations. 

We have compared our measurements with the theoretical model derived in [60] from the multivariate 
halo bias expansion of [25], Eq. (1). At large scales, the confirmation of the validity of the tree- level 
approximation for the halo bispectrum obtained from tliis pcrturbative expansion of the halo density 
is one of the main results of this work. We have studied both the halo bispectrum and the matter- 
matter-halo cross bispectrum, where the lower shot-noise allows for a more accurate comparison between 
predictions and measurements. The value of the constant linear and quadratic bias parameters both 
for Gaussian and non-Gaussian initial conditions are fitted to the measured halo power spectra and 
simultaneously to all, large-scale triangular configurations of the bispectra and later compared to their 
theoretical expectations from the peak-background split approach. 

We have found that the model discussed in Section III provides a quite accurate description, at the 
10% level, at large scales, i.e. k < 0.07/iMpc~^, for almost all triangles of any shape, both for the 
halo bispectrum measured in simulations with Gaussian initial conditions and for the correction to the 
halo bispectrum due to local primordial non-Gaussianity. Since, as a first step, we fit for the constant 
bias parameters and their non-Gaussian corrections, such results signifies that the model presents all 
relevant functional dependencies on the triangular configurations necessary to describe the specific shape 
dependence of the halo bispectrum resulting from nonlinearities in the gravitational evolution and in the 
bias relation between the halo and matter distributions and from the peculiar correlations induced by local 
non-Gaussianity. While the large-scale agreement between the simple tree- level, local bias model and 
numerical results for the halo bispectrum is an established fact for Gaussian initial conditions, recently 
confirmed for instance by [72, 76], the agreement of the model with the additional contribution to the 
halo bispectrum due to local non-Gaussianity, is, on the other hand, not trivial. In fact, the model 
describes such contribution alone by means of up to eight distinct terms, see Eq. (47), each characterized 
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by different scale and shape dependences. We show that for generic triangles none of these terms can 
provide on its own an accurate description of the non-Gaussian effects on the halo bispectruni which is 
rather given by the sum of several different contributions, and, for negative values of quadratic bias, also 
by relevant cancellations between them. This is true, in particular, for squeezed triangular configurations, 
where most of the signal from local models of PNG is concentrated. Interestingly, for such triangles the 
validity of model can be extended to the mildly nonlinear regime, 0.1/iMpc"^ < fc < 0.2/iMpc"^ 

The specific choice adopted for fitting procedure of the bias parameters (where the linear bias parame- 
ters are determined from power spectrum measurements) does not allow for a great freedom to adapt to 
the data, proving to a large extent the predictivity of the model. We compare as well the best-fit value for 
the linear and quadratic bias parameters and their non-Gaussian, scale-independent corrections, to their 
predictions in the context of the peak-background split approach finding a broad, qualitative agreements 
in their mass dependence. 

Finally, we perform a Fisher matrix analysis to compare the ability of power spectrum and bispectrum 
measurements to constrain the non-Gaussian parameter /jvl, and, more importantly, to quantify the 
possibilities given by their combined analysis, a necessary step toward a full exploitation of the data 
available in future large-scale structure surveys. Under the strong assumptions of neglecting the effects 
of covariance and window functions, we show that a combined power spectrum and bispectrum analysis 
can improve over the power spectrum alone by a factor of a few for a very large- volume redshift survey. 
At the same time, the bispectrum can provide a fundamental confirmation of any power spectrum result 
leading, in perspective, to the large-scale structure as a robust test for the initial conditions, with expected 
constraints of the order of those achievable by GMB observations, that is A/jvl ~ few. 

For this to be possible a substantial amount of work is still in order. This paper represents the first 
detailed study and test of a viable model for the halo bispectrum at large-scale. To achieve the accuracy 
necessary to place large-scale structure and CMB observations on the same footing with respect to their 
ability to constrain non-Gaussian initial conditions further investigations are required. In the first place, 
the discrepancies observed between the values of the bias parameters determined from power spectrum 
and bispectrum measurements in [72, 73] might hint at additional contribution, maybe due to nonlocal 
effects, relevant even for Gaussian initial conditions: in general a more accurate determination of the 
properties of halo and galaxy bias will be needed. In the second place, only a comprehensive study of the 
covariance properties of the halo bispectrum in combination with selection function effects can provide 
robust forecasts for the constraints on primordial non-Gaussianity expected from galaxy bispectrum 
observations. Our work will hopefully provide strong motivations for future investigations along these 
directions. 
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Appendix A: Evaluation of matter correlators in Perturbation Theory 

The matter correlators appearing in the halo power spectrum and bispectrum expressions are computed, 
respectively, at 4th and at 6th order in the linear matter density in Eulerian Perturbation Theory 
(EPT) . We refer to the review [77] and references therein for an introduction to cosmological perturbation 
theory. More recent reviews focusing on non-Gaussian initial conditions and higher-order correlators of 
the Large-Scale Structure can be found in [50, 78]. We notice that promising resummation approaches 
in EPT such as those of [79-81] can be extended to non-Gaussian initial conditions [82, 83], leading to 
more accurate predictions than those considered here. 

In addition to the standard loop-corrections to the matter power spectrum and bispectrum wc explain 
here in detail the corresponding nonlinear correction for the cross-power spectrum and bispectrum Pssg 
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and BssSa involving both the nonlinear and non-Gaussian density contrast 5 and its linear and Gaussian 
counterpart Sq. We adopt the notation of [58] for the EPT contributions. 

In Fourier space, the perturbative solution for the nonlinear matter overdensity 6^ is expressed by the 
series [77] 

<5k = 4'^+<5f +4^"^ + ..., (Al) 
where S^^^ = 6o is the linear solution (here non-Gaussian) and 

4")^ Jd%...d^qMciu..-,cin)5^^...S^}, (A2) 

where i^„(qi, q„) is the symmetrized kernel for the n-th order solution. From this expansion one can 
derived in turns perturbative solutions for matter correlators, once the initial conditions, i.e. the initial 
correlators are specified. 



1. Matter power spectra 

Up to 4-th order in Sq and at linear level in /atl, the nonlinear matter power spectrum Pg defined as 
{Si.M = SD{ki2)P6 is given by [18, 84, 85] 

Pg = Pii + Pi2 + P22 + Pl3 + 0{Sl f%L), (A3) 

where, Pn = Pq is the linear matter power spectrum, while the other terms correspond to one-loop 
corrections given by 

Pi2(fc) = 2 yAf2(q,k-q) Po(fc,9,|k-q|), (A4) 

Pi3(fc) = 6Po(fc) y"dV3(k,q,-q) Po(g), (A5) 

P22(fc) = 2 Jd^qFi{ci,k-ci) Po{q) Po(|k - q|). (A6) 
The Gaussian component is 

Ps,G = Pll+P22 + Pl3 + 0{5l), (A7) 

while the non-Gaussian correction is simply 

AAPs,NG = P12 + 0{6l, flL)- (A8) 

The cross-power spectrum Pgg^ is defined as (i^ki^o.ks) = ^d(M-12)P5So so the perturbative expansion 
applies only to one field in the expectation value. We denote the perturbative contributions separating 
the order of correction in 5 from the linear term so that ((^^^j^^o.ka) = <5D(ki2)Pi,i. The EPT expansion 
for Pgs^ is therefore given by 

Pss, = Pi,i + P2,i + P3,i + 0{6lf^L), (A9) 
where Pi^i = Pn = Pq is the linear matter power spectrum, while the 1-loop corrections are given by 

P2,i(fc) = yrfV2(q,k-q) So(fc,(Z,|k-q|) = ^Pi2, (AlO) 

i^3,i(fc) = 3 Po(fc) yd3gP3(k,q,-q) Po(g) = ^Pis- (All) 
The Gaussian component is simply 

P66o,G = Pl,l+P3,l+0{Sl). (A12) 
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while the non-Gaussian correction is APssa,NG = P2,i + ^{^o^ Jnl) which is neglected since it enters at 
second order in in the halo power spectra. 

The evaluation of the matter-matter-halo bispectrum Bgsh and of the halo bispcctrum Bi^, involves 
products of power spectra that, for consistency with the matter bispectra calculation, need to be calcu- 
lated up to 6th-order in the linear density field. In particular 

Ps,G{kl)Ps,G{k2) = Pll{kl)Pll{k2) + [Pll(fci)P22(fc2) + P22(fcl)Pll(fc2)] + 

[Pii(fci)Pi3(fc2) + -Pi3(fci)Pii(A;2)] + 0{5l) , (A13) 

and 

APs,NG{ki)PGik2) = Pi2{ki)Pn{k2) + 0{dl, f^^) • (A14) 
Finally, for the products involving Pggg we have 

PsSo,G{kl)Ps,Gik2) = Pl,l{kl)Pn{k2) + Pl,l{h)P22{k2) + 

[Pi,i(fci)Pi3(fc2) + P3Aki)Pii{k2)] + 0{5l) . (A15) 
2. Matter bispectra 

We evaluate all bispectra (and products of power spectra) up to 6-th order in the linear density field, 
^o- For non-Gaussian initial conditions this implies [58, 86], 

Bs = Bill + Bli2 + ^122 + ^122 + -Bus + -Bus + -B222 + -B123 + -^123 + -^114 + ^Unl^ ^o)>(A16) 

where Bm = Bq is the initial bispectrum and 

= 2 F2(ki,k2) Po{ki) Po(fe) + 2 perm., (A17) 

is the other tree-level contribution, while the 1-loop corrections are given by 

5^22 = 2Po(fci) i^2(ki,k3) Jd^qF2{<i,k3-ci) Bo{k3,q,\k3-q\) + {k3^k2) + 2 perm. 

= F2(ki,k2)[Po(fci) Pi2(fc2) + Po(A:2) A2(fci)]+2pcrm., (A18) 
^(22 = 4 jd^'q P2(q,k2-q) P2(ki+q, kz-q) ^0(^1, <?, |ki+q|) Po(|k2-q|) + 2 perm., (A19) 

Bfi3 = 3Bo{ki,k2,k3) yrf3g^3(k3^q,-q)Po((Z) + 2perm., (A20) 

Biis = 3Po(A:i)y"AP3(ki,q,k2-q)Po(fc2,(Z,|k2-q|) + (fci^fc2) + 2perm., (A21) 

SI22 = 8y"rf35F2(-q,q+ki)P2(-q-ki,q-k2)F2(k2-q,q)Po(9)Po(|ki+q|)Po(|k2-q|), (A22) 

= 6Po(/ci)y"rf3gP3(ki,k2-q,q) P2(k2-q,q)Po(|k2-q|) Po(g) + 5 perm., (A23) 

BiL = 6 Po(A;i) Po(fc2) P2(ki,k2) J d^q F3(ki,q,^) Po(5) + 5 perm. 

= P2(ki, k2) [Po(fci) Pi3(fc2) + Po(fc2) Pi3(fci)] + 2 perm., (A24) 
Pi\4 = 12Po(fci) Po(fc2)yrf^9P4(q -q -ki,-k2) Po(g) + 2 perm.. (A25) 

Specifically, the one-loop contributions present because of non-Gaussian initial conditions are all the 
fifth-order terms Pi22> -^1221 ^113 -Bii3> which depend on the initial bispectrum Bq. The Gaussian 
component to the matter bispectrum is therefore given by 

Bs,G = Bi,2 + 5^2 + BI2S + Pi'l3 + Sn4 + ^^(^o), (A26) 
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while the non-Gaussian correction is 

A-B^^jvG = -Bill + -B122 + -^122 + -Bus + -Bus + 0{5l, f%i), 



(A27) 



Similarly to the cross power spectrum P^^p, for the cross-bispectrum B^sso, defined as (i^ki^ka'^o.ka) = 
Sd{^123)BssSo wc denote the perturbativc contributions separating the order of correction in S from the 
linear term so that {^^l^kz ^0,^.3} = <5D(ki23)-Bij,i- Here the comma indicates that permutations with the 
third index are not considered. We have 

-85550 = -Bll,! + ijfa,! + -^22,1 + -^22,1 + -8^3,1 + -B(|,i + +-828,1 + -823,1 + -8l4,l + ^'(/^i , (Jg XA28) 

where Bn^i = Bm = Bq is the initial bispectrum and 

5(2,1 = 2 F2(ki,k3) Po(fci) Po(fc3) + 2 P2(k2,k3) Po(fc2) Poih), (A29) 
is the other tree-level contribution, while the 1-loop corrections are given by 

^22,1 = Po(fc3)[P2(kl,k3)Pi2(fci)+P2(k2,k3)Pi2(fc2)], (A30) 

4 fd^'q P2(q,k2-q) P2(k3+q, k2-q)So(fc3, |k3+q|,«) Po(|k2-q|), (A31) 



-^22,1 



B' = 3Bo(fci,fc2,fc3) 



-813,1 



Jd^'q P3(ki,q,-q)Po(g) + Jd^'q P3(k2, q, -q)Po(g) 
'6Po{h)jd\iF^ (ki , -q, kg +q)Po (^3 , [kg +q| ) + 
iPo{ki)jd'^q P3(k3,-q,ki+q)Po(fci,?,|ki+q|) + (fci f^fc2), 
-B|3,i = 6 Poik3)Jd^q P3(k3, ki -q, q) P2(ki-q, q) Po(|ki-q|) Po{q) + {h O ^2 
Bill = 6 Po(fci) Po(fc3) P2(ki,k3) J d^'q P3(k3,q,^) Po(g) + (fci ^ ^2) 



= Pl3(fc3) [P2(kl,k3)Po(fcl) +P2(k2,k3)Po(fc2)] : 



(A32) 

(A33) 

(A34) 

(A35) 



B' = 12Po(fc3) 



Po(fci)y'A-F4(q -q -ki ,-k3) Po(g) + Po(fc2)y'A-p4(q -q -k2 -k3) Po(g)(A36) 



In our approximation, we are only interested in the Gaussian component of -85550 given by 

B555o,G = -Bl2,l + -823,1 + -823,1 + -8u,l + ^{JnL^ ^o)- 



(A37) 



We finally notice that the evaluation of the bispectrum contributions for a given triangle takes into account 
in part the effect of the finite bin Ak defining the wavenumbers fci, k2 and ^3. A detailed explanation of 
the procedure is given in Section 3.1 of [61]. 
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